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FOREWORD 


Auburn Research Foundation submitted a proposal which resulted 
in Contract NAS8-28262 being awarded on March 16, 1972. The contract 
was awarded to the Auburn University Engineering Experiment Station 
by the George C. Marshall Space Flight Center, National Aeronautics 
& Space Administration, Huntsville, Alabama, and was active until 
June 15, 1973. 

This report is the final technical report of the work accomplished 
by the Electrical Engineering Department, Auburn University, in the 
performance of the contract. 
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SUMMARY 


This report contains the results of research into the effects on 
system operation of signal quantization in a digital control system. 

The investigation considered digital controllers (filters) operating 
in floating-point arithmetic in either open-loop or closed-loop systems. 

An error analysis technique is developed, and is implemented by a digital 
computer program that is based on a digital simulation of the system. 

As an output the program gives the programing form required for minimum 
system quantization errors (either maximum or rms errors), and the maximum 
and rms errors that appear in the system output for a given bit config- 
uration. The program can be integrated into existing digital simulations 
of a system. 
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I. INTRODUCTION 


The introduction of a digital controller (filter) into a continuous 
data system presents problems to the design engineer that do not exist 
with the use of analog controllers. A major problem in the design of 
digital control systems is the determination of the effects, on system 
performance, of signal quantization within the digital controller. This 
report presents the results of an investigation into the determination of 
the quantization errors, for filters using floating-point arithmetic, 
and the development of design techniques to minimize these errors. 
Throughout this report the terms digital filter and digital controller 
will be used interchangeably. 

A problem in the implementation of a digital filter is the choice 
of the programing form (method of programing) used to realize the filter. 
Generally the use of different programing forms leads to different system 
error magnitudes, caused by signal quantization within the digital filter. 
A considerable amount of research has been published on this topic [See 
References and Bibliography]. In [2], a technique was reported for 
choosing programing forms for filters using fixed-point arithmetic. 

This report presents a technique for choosing programing forms for 
filters using floating-point arithmetic. The research listed in the 
References and Bibliography is concerned generally with digital filters 
in an open-loop configuration. The error analysis techniques require 
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the use of transfer functions, which is not a problem for low-order 
open-loop filters. However, for high-order digital control systems, 
the development of the required transfer functions can be a major 
undertaking. The technique developed in this report is based on a 
digital simulation of the closed-loop system, and thus the pulse 
transfer function of the continuous parts of the system are not required. 
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II. FLOATING-POINT ARITHMETIC 

In this chapter quantization errors that result from arithmetic 
performed in a floating-point format in digital devices are investigated. 

In Chapter III, the results will be applied to the analysis of digital 
control systems to determine system errors resulting from the quantization. 

Floating-Point Format 

In floating-point arithmetic [1], a number ^ is represented as 
the product of two terms, 

x m = E-F (2-1) 

where a part of the bit configuration of the computer word is used to 
represent E, and the remainder to represent F. The term E is the 
exponent and is of the form 2^ for a base 2 computer, 16^ for a base 
16 computer, etc. , where y is a signed integer. The bit configuration 
for E yields the value of y. The term F is the fraction, and is normally 
set such that 

1/2 < F < 1, (2-2) 

for a base 2 computer, if the base of the computer is 2 k , 
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(2-3) 


l/2 k < F < 1 

The number zero is a special representation. The bit configuration 
for F yields F directly, where the first bit in F represents 1/2, the 
second bit 1/4, etc. 

Let s be the number of bits assigned to the exponent, and t the 
number assigned to the fraction (excluding the sign bit for the 
fraction). Also let 

s + t = n 


Thus there is a total of n 4- 1 bits in the computer word configuration, 
with the additional bit used to give the sign of the number represented. 
The maximum magnitude of the exponent for a base 2 computer is 


E 

m 



( 8 - 1 ) 


- 1 ] 


(2-4) 


and, for a base 2 computer, is 


kf 2 1 

E m = 2 [2 l] 


(2-5) 


The factor (s-1) appears since one bit of s must be used to give the 
sign of the exponent. The maximum magnitude of the fraction, F, is 
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( 2 - 6 ) 


F = 1 - l/2 t 
m 

Thus the maximum magnitude that can be represented by the n bits is 

= ? k[2(S ” 1) - 11 [l-l/2 t l « 2 kf2 (s_1) - 1] (2-7) 

for a base 2 k computer. Table 2-1 lists the numbers that can be 
represented in a base 2 computer by a bit configuration with n equal 
to five. In this configuration, s is equal to three, and is the first 
three bits. Then t is two, and F is represented by the last two bits. 

If F satisfies (2-2), the fourth bit from the left in the configuration 
is always 1. These values are indicated by an asterisk. The bit 
representation for zero is also shown by an asterisk. The truncation 
quantization characteristic is shown in Figure 2-1 for this case. 

Quantization Errors 

The characteristics of the quantization errors will be determined 
in this section [3], Let Xjjj be the floating-point machine representation 
of x. Then 

Xm = Q f £ [x] = 2 k Y. F , (2-8) 

where Q^t*] indicates the floating-point representation 
(which is quantized) of the number. Suppose that truncation, as shown 
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TABLE 2-1 

EXAMPLE OF QUANTIZATION 


bit 

configuration 


floating-point 

number 


11111 

11110 

11101 

11100 

11011 

11010 

11001 

11000 

10111 

10110 

10101 

10100 

10011 

10010 

10001 

10000 

01111 

OHIO 

01101 

01100 

01011 

01010 

01001 

01000 

00111 

00110 

00101 

00100 

00011 

00010 

00001 

00000 


6 * 

4* 

2 

0 

3* 

2 * 

1 

0 

3/2* 

2 / 2 * 

1/2 

0 

3 / 4 * 

2/4* 

1/4 

0 

3/32* 

2/32* 

1/32 

0 

3/16* 

2/16* 

1/16 

0 

3/8* 

2 / 8 * 

1/8 

0 

3/4* 

2/4* 

1/4 

0 * 
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in Figure 2-1, is used in quantizing x. Then the maximum magnitude of 
the quantization error is seen to be 


e m = 2^11/2^ 


(2-9) 


and this error is always negative. For roundoff quantization, the 
maximum error is one-half that of (2-9). Then, from (2-8) and (2-9), 


e m " x m ' 2 


-t 


/F 


( 2 - 10 ) 


Thus the quantization error is maximum if F is minimum. For a base 
2 computer, from (2-2) , 


e m = x m' 2 “ (t ' 1) (2-11) 

For a base 16 computer, 

e m = (2-12) 

lc 

Thus for a base 2 computer, where 

l/2 k <_ F < 1, (2-13) 


then 



e m = V 


-(t-k) 


(2-14) 


*2 


For round-off quantization, the maximum error is one-half that given 
in (2-14). It is necessary that x and Xjjj be approximately equal, or 
else all calculations are meaningless. Thus x ra may be replaced with 
x in (2-14) , with very little resulting error. Thus 

e m = x-2‘ (t ' k) (2-15) 

In general, the error introduced by truncation quantization is 

e = 2^6 ; 8 < 1/2* (2-16) 

as is seen from Figure 2.1. Thus, from (2-8), 

e = x^g/F = 6 (2-17) 


where 


S 


< 


2 -(t-k) 


(2-18) 


for truncation, and 


( t-k+1) < 6 < 2 -(t-k+l) 


(2-19) 
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for round-off. Generally in this type of analysis, it is assummed that 
all errors in the ranges given by (2-18) and (2-19) are equally likely 
to occur. This approach was first taken in this investigation, but the 
results were grossly in error when compared to actual errors obtained 
from simulation. Thus it was found to be necessary to use a more accurate 
representation. For floating-point quantization, the probability density 
functions for round-off and for truncation are given in Figure 2-2 [3] . 
These density functions will be used in the analysis developed in 
Chapter III. 

Signal quantization errors occur during four operations in a digital 
filter. These four operations are: (1) analog- to- digital conversion, 

(2) digital- to- analog conversion, (3) multiplication, and (4) addition. 
The errors in (1) and (2) are described above. For multiplication and 
addition [1], 

Q f& IxyJ = (xy)(l + 6), (2-20) 

and 

Q fJl Ix + y] = (x + y)(l + «), (2-21) 

where 6 is given by (2-18) and (2-19). However, care must be exercised 

[3] in the utilization of (2-21). Consider the addition of two numbers 
in base 10, with the numbers limited to three decimal places. 
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.115 

X 

10 3 

.115 

X 

10 2 

.173 

X 

10 3 

.173 

X 

<n 

o 

*H 

.288 

X 

10 3 

00 

H 

5 x 

10 3 


In the first case, with the exponents equal, these is no error in the 
addition. In the second case, with the exponents unequal, there is 
a quantization error. Thus, in a digital filter, if the exponents of 
the numbers being added are equal, there will probably be no error in 
the addition. If the exponents are unequal, the density functions for 
the errors are given in Figure 2-2. 



(a) Truncation ( b > Round-off 

Figure 2-2. Probability density function for 6. 
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III. SYSTEM ERROR ANALYSIS 


In this chapter, the error models derived in Chapter II will be 
utilized to develop a technique for determining system errors due to 
quantization in a digital control system. The technique will be 
implemented using a computer simulation of the control system. 

System Errors 

It was shown in Chapter II that the quantized representation of a 
number in floating-point format can be written as 

Qf^ x) = x(l + <5) (3-1) 


or 


x - Q f£ (x) - x6 


(3—2) 


Thus the quantization operation may be modeled as shown in Figure 3-1. 



Figure 3-1. Model of floating-point quantization. 
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Consider now an otherwise linear discrete system that contains a 
single quantizer. Let the 



Figure 3-2. Discrete system containing quantization. 


system be represented as shown in Figure 3-2. The system (b) of 
Figure 3-2 can be considered to be linear if the sequence x(k)6 k , 
the error at the k th sampling instant, is determined in advance. The 
output y(k) then can be expressed as 

y(k) = y r (k) + y q (k) , (3-3) 

where y r (k) is the response from the input r(k), and y q (k) is the 
response from x(k)6 k . Thus y q (k) is the error in the output due to the 
quantization. Let G(z) be the transfer function from r(k) to x(k) . 

Then 

X(z) = R(z)G(z) = x(0) + x(l)z” 1 + ... (3-4) 
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Let 


X 6 (z) = 5 q x( 0) + 5jX(l)z 1 + 6 2 x(2)z -2 + .... (3-5) 

Also let H(z) be the transfer function from the error source x(k)S k to 
the output y(k). Then 

Y q (z) = H(z) X (z) 

= 6 0 x(0)H(z) + 6 1 x(l)z" 1 H(z) + 6 2 x(2) z'hiiz) + ... (3-6) 

At the k th sampling instant, 

y q (k) = 6 0 x(0)h(k) 4- 6ix(l)h(k-l) + ... + 6 k x(k)h(0) , (3-7) 

where (h(k) } is the impulse response from the error source to the output, 
given by 

H(z) = h(0) + h(l)z’ 1 + h(2)z" 2 + ... (3-8) 

Assume that the quantization is round-off. Then 6^ is given by (2-19), 
which is repeated as (3-9) below. 

_ 2 -(t-k+!) < 6 2 -(t-k+l) (3-9) 
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Thus it is seen from (3-7) that the maximum magnitude of the error in the 
output is given by 


yqmaxOOl = 6 max [ |x(0)h(k) | + |x(l)h(k-l) | + ... + |x(k)h(0)|], 

(3-10) 


where 


1 1 

°max 


2 -(t-k+i) 


For a base 2 machine. 


(3-11) 


«max - 2-t (3-12) 

Consider now truncation quantization. For truncation, the error 
is alway negative or zero. Thus, in (3-7), each 6^ is either negative 
or zero. The maximum possible value for yq(k) is obtained from (3-7) 
by allowing 6^ t0 assume its maximum magnitude, and first summing all 
positive terms, and then summing all negative terms. The maximum 
possible magnitude of y q (k) is then the larger of the two sums, in 
magnitude. It is seen that larger of the two sums is at least one-half 
the sum of (3-10). However, it is to be recalled that <5 max for trunca- 
tion is twice that for round-off, for a given t and k. 
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To determine the mean-square error, consider equation (3-7). The 
output y q (k) can be considered to be the output of a system whose input 
is given by 

A (z) = <5 q + 5^z ^ + <$ 2 2 ^ + ... > (3—13) 

and whose tranfer function is given by 

F(z) = h(0)x(k) + h(l)x(k-l) z” 1 + h(2)x(k-2)z“ 2 + ... 

- f(0) + f(l)z" 1 + f(2)z“ 2 4- ... (3-14) 

The 6 ± of (3-13) are assumed to be independent, and to have the density 
functions given in Figure 2-2. The following development will be for 
both round-off and truncation quantization. Let m^ be the expected 
value of and m 2 be the second moment, i.e. , 

E£6 ± ] = m i; E[6^] = m 2 ; for 1 (3-15) 

Consider now (3-7), (3-13), and (3-14). The expected value of the output 
error resulting from a single quantizer is 


E[y q (k)] = ElS 0 f(0) + 6if(l) + ... + S k f(k)] 
k 

= mj. 1 (3-16) 

i=0 



and m^ 


are derived in Appendix A. 
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The mean-square output error is given by 


E[yq(k)] = E[{« o f(0) + 6^(1) + ... + « k f(k)} 2 ] (3-17) 


Thus 


E[y 2 (k)] = E [^q] f 2 (0) + ... + E[6^]f 2 (k) + 2f(0)f (1) E [6 (0) ]E £6 (!) ] 


+ 2f (0)f (2)E [6 (0) ]E [6 (2) ] + ... + 2f (k-l)f (k)E[6(k-l) ]E[6 (k) ] 

(3-18) 


Or 


k k-1 k 

E[y 2 (k)] = m 2 I f 2 (i) + 2m f J £ f(i)f(j) (3-19) 

H i=0 i=0 j=i+l 


But 


k-1 


2 l 

i=0 


k 

I f(i)f(j) = [ 
j=i+l 


k 9 

l f(i) ] 

i=l 


i f 2 a) 

i=l 


(3-20) 


Then (3-19) becomes 


E[yq(k)] 


k k 

l f 2 (k)[m 2 - mf] + I l f(i) ] 2 m 2 
i=0 i=l 


(3-21) 


Consider now the case that n quantization points contribute to the output 
error. Let y t ( ^9 be the total error in the output, and let yq^(k) be that 
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part of ytOO due to the i 1 -* 1 quantization point. Then, the mean output 
error is 

Ely t (k)] = E[y ql (k) + y q2 (k) + ... + y qn (k)] 

= Ely ql (k)] + Ely q2 (k)] + ... + E[y qn (k)], (3-22) 

where E[y^(k)] is given by (3-16). The mean-square output error is 
given by 

Ely^(k)] = El{y ql (k) + ... + y qn (k)} 2 ] - EHy ql (k)} 2 ] + ... 

+ E[{y qn (k)}2] + 2Ely ql (k)]E[y q2 (k)] + ... 

+ 2 E[y q(n -i)(k)]E[y qn (k)J (3-23) 

The terms of (3-23) are given by either (3-21) or (3-16) . 

Error Calculations by Digital Simulation 
Both the maximum error of (3-10) and the mean-square error of (3-21) 
can be obtained from digital simulations of the system. First, with the 
input of the system of Figure 3-2 equal to the desired input {r(k)} and 
all initial condition set of zero, the sequence (x(k)} can be calculated 
and stored. Next, with zero input and zero initial conditions, the 
impulse response sequence {h(k)} can be obtained by applying are input 
6qx(0) = 1, and 6jx(i) = 0, i > 0. The system output sequence is then 
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the desired impulse response sequence {h(k)>. The indicated sums of 
(3-10) and (3-21) can then be evaluated, yielding the maximum error 
and the mean-square error for the system. The total maximum error 
is then the sum of the maxi mum errors from each quantizer, and the total 
mean-square error is given by (3-23) . 

Computer Evaluation of Errors 

A computer program has been developed for determining system 
quantization errors in a digital control system, with the digital filter 
realized by various programing forms. The utilization of this program 
will allow the choice of the programing form that will yield the smallest 
system error due to quantization. This computer program will be 
described in this section. 

All programing forms in the computer program realize the transfer 
function 

2 

a^z + ai z + a n 

D(z) = 0 

z^ + b^z + bQ 

Consider first the canonical programing form, shown in Figure 3-3. It 
is assumed that this filter form is connected in a digital control system. 
Each point at which quantization occurs is indicated by an The input 

and output quantization points are omitted, and are considered separately, 
since the errors from these points are the same for all programming 
forms. Let X2(z) be the signal at quantization point 1 as indicated in 
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Figure 3-3. Further, let H^(z) be the transfer function from the 
error source at quantization point 1 to the system output. Then the 
errors resulting from quantization point 1 are obtained using X2(z) 
and H^(z) in (3-10) and (3-21). The errors resulting from quantization 
point 2 are obtained using biz _1 X2(z) and -H 1 (z) in (3-10) and (3-21). 

All necessary signals and transfer functions to determine system errors 
from the eight points of quantization are given in Figure 3-3. 

The sums indicated in (3-10) and (3-21) are evaluated using the 
computer program given in Appendix B of this report. The value of k 
in these equations is N1 of the program, and must be given in the main 
program. The canonical programing form is simulated in subroutine 
F1L2. Three different simulations of the system are required. In the 
first simulation (JJ = 1), the sequence {x2(k)> of Figure 3-3 is obtained 
and stored. Note that the system input is R, and is a unit step in 
this case. In the second simulation (JJ = 2), the impulse response 
{hj_(k)} is obtained, and £f(i) and Sf^(i) are evaluated for quantization 
points 1, 2, 3, and 4. In the third simulation (JJ = 3), the impulse 
response is obtained, and the required sums are evaluated for 

quantization points 5,6,7and 8 

The quantization error at any point is assumed to be statistically 
independent of that at any other point. Thus the maximum system error 
for the filter form is the sum of the maximum system errors for each 
quantization point, and the mean-square system error is given by (3-23). 
The total system errors are evaluated using these relationships. The 
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signal at 

transfer function 

error point 

error point 

to output 

3 

X2 (z) 

H x (z) 

4 

b 1 z' 1 X2(z) 

-^(z) 

5 

b 0 z' 2 X2(z) 

-^(z) 

6 

(b 1 z“ 1 +b 0 z" 2 )X2(z) 

-H (z) 

7 

a 2 X2(z) 

H 2 (z) 

8 

a^z”^X2 (z) 

H 2 (z) 

9 

a Q z“ 2 X2(z) 

H 2 (Z) 

10 

(a^+a^z -1 ) X2 (z) 

h 2 ( Z ) 


Figure 3-3. Canonical form. 





program prints out the maximum error and the root-mean-square error 
for each filter form with each divided by 6 max of (3-10) and (3-19) . 
Since no errors occur in multiplication if a signal is multiplied by 
a coefficient equal to unity, logic is included in the program to zero 
the multiplication errors for this case. To obtain the errors for a 
given system, the results of the computer program must be multiplied 
by <$ 

max 

Consider now the modified canonical programing form. This form 
is shown in Figure 3-4, and is realized as subroutine F1L1 in the com- 
puter program. This form is also used in determining the errors from 
both input and output quantization. For this programing form the 
sequences {X2(k)} and {EI(k)> are calculated and stored for JJ = 1. 

For JJ = 2, the impulse response from X2 to the system output is 
obtained. For JJ = 3, the impulse response from the filter output is 
obtained; and for JJ = 4, the impulse response from the filter input 
is obtained. These responses are used in the manner indicated in 
Figure 3-4 to calculate system errors. 

All filter forms consider are listed in Table 3-1. It is noted 
that the direct form gives the same system errors as does the canonical 
form [4], and thus is not considered separately. The parallel form can 
realize filters containing only real poles, and the XI and XII, forms 
can realize filters containing only complex poles. Logic is included 
in the program to insure that these forms are considered only at the 
appropriate times . 
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error point 

signal at 
error point 

transfer function 
to output 

3 

X2 (z) 

H^z) 

4 

b x z~h2(z) 

-H^z) 

5 

b(f - 2 X2(z) 

-HjCz) 

6 

4 + 5 

-Hx (as) 

7 

a 2 EI(z) 

H 2 (z) 

8 

(a i -b 1 a 2 )z” 1 X2 (z) 

H 2 (z) 

9 

(a 0 -boa 2 ) 2 '” 2 X 2 (z) 

H 2 (z ) 

10 

7+8 

H 2 (z ) 

2 

9 + 10 

H 2 (z ) 

1 

EI(z) 

H 3 G 2 ) 


Figure 3-4. Modified canonical form 
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TABLE 3-1 


FILTER FORMS 


filter form 

Figure 

Subroutine 

modified canonical 

3-4 

FlLl 

canonical 

3-3 

F1L2 

modified direct 

3-5 

F1L3 

direct 

errors same 

as canonical 

modified standard 

3-6 

F1L4 

standard 

3-7 

FIL5 

parallel 

3-8 

F1L6 

XI 

3-9 

F1L7 

XII 

3-10 

F1L8 
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El EO 




signal at 

transfer function 

error point 

error point 

to output 

3 

• X2(») 

%(z) 

4 

b 1 X2(z) 


5 

b 0 z -1 X2.(z) 

-z' 1 H 1 (t ! ) 

6 

4 + 5 

-z‘ 1 H I (z) 

7 

32EI(z) 

H X (Z) 

8 

a x £ _1 EI(z) 

V«> 

9 

-2 

a Q z EI(z) 

H x (z) 

10 

7 + 8 

H^z) 

11 

9 + 10 

%(z) 


Figure 3-5. Modified direct form. 
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E 



3 X2(z) Hjtt) 

4 a 0 EI(z) H 1 (z) 

5 b 0 {z^ 1 [X3(z)+a 1 EI(z)]+a 2 EI(z)} -H^z) 

6 X3(z) H 2 <») 

7 a 1 EI(z) H 2 (») 

8 a 2 EI(z) z."' 1 H 2 ( a ) 

9 z’ 1 [X3(z) + a 1 EI(z)J + a 2 EI(z) z~ 1 H 2 (z) 

10 b 1 {z~ 1 [X3(z) + a-jEKz)] + a 2 EI(z)} ” H 2 (Z) 

11 6+7 H 2 (2) 


Figure 3-6. Modified standard form. 
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3 

X2(z) 

H x (z) 

4 

b l Z ~ lx 2 (Z) 

-H 1 (z) 

5 

^ a 0‘" a 2 b 0 _b l (a l“ a 2 b l )]EI(a) 

H x (z) 

6 

b 0 [(a - a 2 b i )r 1 EI(z) + z" 2 X 2 (z)] 

-%(*) 

7 

4 + 6 

- H l(z) 

8 

(a 1 -a 2 b 1 )EI( z ) 
-1 , „ 

H 2 (z) 

9 

z X2(z) + (a 1 -a 2 b 1 )EI(z) 

H 2 (z) 

10 

a 2 EI(z) 

H 3 (z) 


Figure 3-7. Standard form. 
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EO 

-O 


error point 


signal at 
error point 


transfer function 
to output 


3 

4 

5 

6 

7 

8 
9 

10 


gjEKz) 

XI (t) 

p 1 z" 1 Xl(z) 

g 2 EI(z) 

X3(z) 

-1 

$2 Z x3 (*> 
t^XKz) + z _1 X3(z) 

a 2 EI(z) 


H x (z) 

H x (z) 

-%(*) 

H 2 (z) 

h 2 (z) 

-H 2 (2) 

V z) 

H 3 (z) 


Figure 3-8. Parallel fora. 






3 

Xl(z) 

-1 

H X (Z) 

4 

P 3 z Zl(z) 

-^(z) 

5 

2g 3 EI(z) 

H^z) 

6 

V^z~ l X3(*) 


7 

5-6 


8 

X3(*) 

-1 

H (Z) 
2 

9 

P 3 z X3(z) 

" H 2 (Z) 

10 

2g 4 EI<z) 

h 2 ( z ) 

11 

P 4 z” 1 Xl(z) 

h 2 ( z ) 

12 

10 + 11 

H 2 (Z) 

13 

a 2 EI(z) 

H 3 ( z ) 


Figure 3-9. XI form. 
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error point 


signal at 
error point 


transfer function 
to output 


3 Xl(z) H x (z) 

4 P 3 r 1 Xl(z) -H 1 (Z) 

5 EI(z) + P 4 z -1 X3(z) H-^z) 

6 P 4 Z _1 X3(Z) H x (z) 

7 X3(fc) H 2 (z) 

8 P 3 'z _1 X3(Z) -H 2 (z) 

9 p 4 z' 1 Xl(z) H 2 (z) 

10 2g 3 z -1 X3(z) H 3 (Z) 

11 2g 4 Z _1 Xl(z) H 3 (Z) 

12 10 + 11 H 3 (z) 

13 a 2 El(z) H (Z) 


Figure 3-10. XII form 





Higher Order Filters 

In the previous sections, it was assumed that the digital filter 
was of at most second order. In this section it will be shown that 
the program given in the appendix is applicable to higher order filters. 

Consider, for example, a fourth-order filter, with all poles 
complex. This filter should be realized as two second-order sections, 
as shown in Figure 3-11, because of coefficient sensitivities. There 
are two questions to be answered with respect to this filter. First, 
which section should be placed first, and second, which programing form 
should be used for each section? To answer these questions, four 
different computer runs must be made. The filter section D^Cz) must 

2nd order 2nd order 



Figure 3-11. Digital control system. 

be simulated in the program first as shown in Figure 3-11, and then with 
its position reversed with that of D 2 (z). Then the same simulations 
must be made with D 2 (z). An examination of the results of these runs 
will indicate not only the placement of D^(z) and D 2 (z) , but also the 
filter forms required. For the system of Figure 3-11, an addition 
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subroutine must be added to the program to simulate the additional 
filter section. 


Example 

As an example, consider the system of 


, , glaat 


V 

> 

z 2 -l . 25z +.32 

^ 

z + .5 ! 


— > 

\ 

/ * 



z^ -1.6z + .63 










Figure 3-12. System for example 


Figure 3-12. This system is simulated in the program of Appendix B, 
and the predicted errors from quantization were calculated. To check 
these results, the system was simulated with first the filter in single 
precision and the remainder of the system in extended precision and then 
with the entire system in extended precision on an IBM 360/50 computer. 

The difference in the outputs of these two simulations is then the system 
quantization error. The IBM 360/50 computer is a base 16 machine, has 
24 bits in the fraction, and used tuncation quantization. Thus, in (2-18), 
t is equal to 24 and k is equal to 4. The maximum error at the point of 
quantization is 

6 m = 2 0.95 x 10" 6 (3-24) 
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Consider first the modified canonical form. The predicted rms 
error for the system was found to be 0.80 x 10 **, with k of (3-21) 
equal to 100 and the input a step function. 1 Sixteen different simulations 
with step-function inputs yielded a total rms error of 0.33 x 10 , 

which is lower than the predicted result. However, it is to be recalled 
that probably no error will occur in the addition of two numbers if 
the exponents of the two numbers are equal. Since the IBM 360/50 is 
a base 16 computer, two numbers can be different by a factor of almost 
16 and still have the same exponent. To approximately compensate for 
this effect, the errors from all additions were set to zero, and the 
predicted rms error was calculated to be 0.28 x 10“^. It is then seen 
that additions do contribute very little to the total system errors. 

Next the canonical form was used in the simulations, with the 
resultant rms error of 0.20 x 10~ 6 . With addition errors included, the 
predicted rms error was 0.27 x 10~ 6 ; and without additional errors, the 

predicted rms error was 0.24 x 10 ”**. 

With the standard form, the result of the simulations was an rms 
error of 0.135 x 10” 6 . With addition errors included, the predicted 
rms error was 0.21 x 10“ 6 ; and without addition errors, the predicted 
rms error was 0.125 x 10 

From the above results it is seen that the developed technique 
yields accurate results. However, problems do occur with the errors 

1 Computer run time for the program of appendix B was approximately 

30 seconds, using the WATFIV version of FORTRAN. 
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caused by addition. In a base 2 computer, these errors would be 
much more likely to occur, since the exponents of the numbers are 
not as likely to be equal. Additional research is needed to develop 
a better technique for the inclusion of errors due to addition. 
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IV Conclusions 


In this report the problem of floating-point quantization errors 
in a digital control system is investigated. A technique is developed 
which yields the maximum possible system output error and the mean- 
square system output error caused by quatization in a digital controller 
operating with floating-point arithmetic. A deterministic system in- 
put is assumed. The technique is implemented by a digital computer 
program, which is based on a simulation of the control system. The 
program requires as input the number base of the digital controller, 
the digital controller coefficients, the system input function, and the 
sampling period for which the errors are desired. The system plant 
must be simulated in a subroutine of the program. As an output, the 
program gives the maximum error and the mean-square error in the system 
output for the digital controller realized by nine different programming 
forms. Thus the programming form can be chosen for the controller that 
yi£lds the smallest errors. 
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APPENDIX A 


In this appendix the first and second moments of the probability 
density functions for both roundoff and truncation errors will be derived. 
Consider first the roundoff density function 



Figure A-l. Roundoff quantization 
given in Figure A-l* Let r = 2 k be the base of the con^uter. The 

first moment is zero. The second moment is 



6 2 f (6 ) d<5 - 



<S 2 f (6)d6 


(A-l) 


Or 


m 


2 



h « dfi + 2 


l 


6 m 


m 


-h r 


(« - 6 m )6 2 d« 


(A-2) 


Since the area under the curve in Figure A-l is equal to unity, then 


h 




(A-3) 
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Substituting (A-3) into (A-2) and evaluating the integrals, we find that 


6 m 2 (r 3 + r 2 + r + 1) 
6r 2 (r + l) 


(A-4) 


Consider now the truncation density function given in Figure A-2. 


The first 


m 


1 




<5f(6)d<5 



- h r 6 

6 m ( r -D ~ 6 m )d6 » 


(A-5) 


where 


2r 

h = (r+l)6 m 


(A-6) 


Evaluation of (A-5) yields 


dm( r2 + r + 1) 

3r(r + 1) 


The second moment is given by 




m. 


= / <S 2 hd5 + f 


m 


-h r 5 


<$ (r-1) (<5 " 
m v 7 


5 m )dd 


(A-7) 


(A-8) 
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Evaluation of (A-8) yields 


<$_ 2 (r 3 + r 2 + r + 1) 
III 2 = Mi l ' 

6r^ (r + 1) 


(A-9) 


Note that (A-4) and (A-9) are the same. However, for roundoff is 
one-half that for truncation. 
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APPENDIX B 


TABLE, b-1 
SYMBOLS IN PROGRAM 


Symbol 

Description 

BK 

base of computer 

BSR 

use to determine if poles of filter 
are complex 

El (I) 

signal amplitude at quantization point 

EE (I) 

£ x(Nl~i)h(i) 

EER 

used in calculating rms errors 

El 

input to filter 

EH 

used in finding impulse response 

EMM 

total maximum error for form 

EM(I) 

£|x(Nl-i)h(i)| 

EO 

output of filter 

ERMS 

total RMS error for form 

ERR 

used to calculate total RMS error for form 
2 

ER(I) 

£[x(Nl-i)h(i)] 

F(I) 

x(Nl-i)h(i) 

10 

10=0 zeros errors from input quantizer 

11 

11=0 zeros errors from output quantizer 

JJ 

number count for simulations for each form 

JJJ 

number of simulations for each form 

JT 

JT=0 for roundoff, JT=1 for truncation 

KK 

used to choose form 

KKl 

initial value of KK 
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Table B-l 
(continued) 


KK2 

MM 

N1 

R 

XM 

XM2 

XM3 

YO 


Final value of KK 

stops program at correct point 
for filter with complex poles 

number of iterations in each 
simulation 

system input 

first moment (expected value) 
for errors 

used to calculate rms errors 
second moment for errors 
system output 
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START 



Figure B-l. Flow-chart for program 
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Figure B-l (continued) 
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3 



IS 

- KK2 


/'ARE \ 

FILTER POLES 
REAL 


Figure B-l (continued) 
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PROGRAM 


DIMENSION E 1 (150) ,62 < 150) ,E3<150) ,FC 131 ,ER( 13),EM< 133 ♦ EE ( 13) 
C0MM0N/C1 / A2 1 A1 , AO, B1 , 80, X0 , XI , X2 ,X3 , X4 , X 5, E I l , E I 2, 6 1 3 
C0MM0N/C2/El*E2,e3,F,ER,EM*EE, II* JJ,KK*.Nl,YO* JT*XM,XM2 
C0MM0N/C3/P 1, P2,P3»P4»GltG2,G3,G4 
COMMON/C4/A22,A11,AOO,B22,B11,BOO,G,X31 ,X4l,X51 

C JT=0 FOR ROUNDOFF, JT=1 FOR TRUNCATION 

JT * 1 

J T = 0 

C BASE OF COMPUTER I S BK 

BK * 16. 

XM IS FIRST MOMENT OF ERROR SOURCE 
XM3 IS SECOND MOMENT OF ERROR SOURCE 
XM2 IS USED IN CALCULATING RMS ERRORS 
XM* ( BK**2*3. *BK+3 • )/(BK**2*3.*BK+2.)/34 

XM3*< BK**3«-4.*BK**2«-6.*BK4-4. )/( BK**3+4. *BK**2+5.*BK+2. 1/6. 
IFUT.NE.l) XM = 0. 

XM2 = XM3 - XM**2 

C COEFFICIENTS OF PLANT 

G = 1.0 
A22 = l. 

All = 0.5 
AOO = 0, 

B22 * -1.6 
Bll = 0.63 
BOO * 0. 

C COEFFICIENTS OF FILTER 

AO * 0.315 
A 1 * -1.25 
A2= 1 • 0 
BO * 0.035 
B 1 = -0.75 
MM * 1 

C TO DETERMINE IF FILTER POLES ARE COMPLEX 

BSR = (Bl**2)/4. - BO 

600 IFIBSR.GT.O.) GO TO 601 

C FORMS FOR COMPLEX POLES 

KK1 = 1 
KK2 = 5 
GO TO 603 

601 CONTINUE 

C FORMS FOR RBAL POLES 

BS SR = SQRTTBSR) 

PI = Bl/2. * BSSR 
P2 = Bl/2. - BSSR 

Gl * ( A 1*P l - A2*P1**2 - AOJ/IPl - P2) 

G2 * (A0-A1*P2 ♦ A2*P2**2)/(P1-P2> 

KK1 = 1 
KK2 = 6 
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GO TO 603 

604 IFIBSR.GE.O. ) GO TO 606 
MM * 2 

C FORMS FOR COMPLEX POLES 
P3 * Bl/2. 

P4 =5 SQRTI-8SRI* 

G4 * A 1 /2 • - P3*A2 

G3 = ( AO - A2*B0 - 2.*G4*P3)/<2.*P4) 

KK 1 = 7 
KK2 * 8 
603 CONTINUE 

C N1 IS TOTAL ITERATIONS FOR SIMULATIONS 

N1 = 100 

C KK USEO TO CHOOSE FILTER FORM 
DO 501 KK-KN1.KK2 

C JJJ SETS NUMBER OF SIMULATIONS REQUIRED FOR EACH FORM 

JJJ = 4 

IF(KIK*EQ«2) JJJ?=3 
IFIKK.EQ.3) J J J = 2 
IF ( KK • EQ« 4 ) J JJ a 3 

C JJ COUNTS THE NUMBER OF SIMULATIONS FOR EACH FORM 

00 50 JJ«1 # JJJ 
YO * 0. 

C ZERO INITIAL CONDITIONS FOR FILTER 

XO = 0. 

Xl-O. 

X2*0, 

X3 a 0. 

X4=0. 

X5 = 0. 

C ZERO INITIAL CONDITIONS FOR PLANT 

X31 = 0. 

X41 » 0. 

X5i = 0. 

C L I 1 » E I 2 » E I 3 USED TO CALCULATE IMPULSE RESPONSES 

Eil = 0. 

£12 a 0 . 

E I 3 = 0. 

C Eim f E2(I)fE3(n ARE SIGNALS AT POINTS IN FILTER 
Eld) = 0. 

EK2) * 0. 

E2( 1) = 0. 

E2I2) = 0. 

E 3 ( 1 ) = 0. 

E 3 ( 2 ) a o. 

JJ*2 YIELDS FIRST IMPULSE RESPONSE 
J J=3 YIELOS SECOND IMPULSE RESPONSE 
J J a 4 YIELDS THIRD IMPULSE RESPONSE 
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I F ( J J . EQ. 2 ) Ell = 1. 

IFUJ.EQ.3) E 12 * l. 

IFIJJ.EQ.4) E 13*1 • 

C JJ=1 YIELDS SIGNAL AMPLITUDES IN FILTER 

IFIJJ.GT.2) GO TO 102 
IFIKK.GT. l) GO TO 10 
ERII) USED TO CALCULATE RMS ERROR 
EE ( I ) USED TO CALCULATE RMS ERROR 
EMI I) USED TO CALCULATE MAXIMUM ERROR 

00 1 1= It 2 

EE ( I ) * 0. 

ERII) * 0. 

1 EMI I) * 0. 

10 DO 2 1*3, 13 

EE { I ) = 0. 

ERII) * 0. 

2 EMI I ) = 0. 

C R IS SYSTEM INPUT 

102 R = 0. 

IF ( JJ.EQ. 1) R =. 1. 

IFIJJ.EQ.2) N 1 * N 1 - 2 

DO 40 11*1, N1 

CALL PL ANT I BI , EO* R ) 

40 CONTINUE 

50 CONTINUE 

501 CONTINUE 

IFIMM.EQ. 1) GO TO 604 
606 CONTINUE 
STOP 
END 
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SUBROUTINE PLANT ( E I « EO,R ) 

C PLANT OF CONTROL SYSTEM 

0 1 MENS I ON El (150) ,E2( 150) ,63(150 > f Ft 13) ,ER( 13), EM( 13>,EE( 13) 
COMMON/C2/El»E2*E3 f F f ER»EM*EEt 1 1 * JJ*KK r Nl » YO f JT , XM , XM2 
COMMON/C4/A22,All,AOO,B22,Bll.BOOtG*X31*X41, X51 
C YO IS PLANT OUTPUT 

YO = A00*X31 * A 1 l*X4 1 ♦ A22*X51 
El * R - YO 
El IS FILTER INPUT 

EO IS FILTER OUTPUT 

IF(KK.EQ.l) CALL FIL1(EI«E0) 

IF ( KK • EQ. 2 ) CALL FlL2{EitEG> 

IF(KK.EQ.3) CALL FIL3(EI«EOt 
IF ( KK *EQ* 4 ) CALL FIL4I6I,E0) 

IFIKK.EQ.5) CALL FIL5(EI,E0) 

I F ( KK « EO* 6 ) CALL FIL6(EI*E0) 

IFIKK.E0.7) CALL F I L7 ( E I * EG) 

IFIKK.E0.8) CALL FIL0<EI,EO) 

X6i = G*EO - B00*X3i - B11*X41 - B22*X5l 
X31 * X41 
X41 * X51 
X51 = X61 
RETURN 
END 


49 



ooo ooo o no noon 


SUBROUTINE FILH£l,EO* 

C MODIFIED CANONICAL FORM 

DIMENSION Ell 150) ,E21 150>,E31 150),F113) ,ER( 13) ,EM(13) ,EE< 13) 
COMMON/C1/A2, A1 , A0,B 1 , BO, XO, XI , X2 , X3 , X4, X5, E I i , E 1 2, E 1 3 
COMMON/C2/El,E2,E3,F,ER,EMiEE, II , JJ, KK,.Ni , YO, JT,XM,XM2 
El IS FILTER INPUT 

Ell USED TO CALCULATE FIRST IMPULSE RESPONSE 
EI2 USED TO CALCULATE SECOND IMPULSE RESPONSE 
E 1 3 USED TO CALCULATE THIRD IMPULSE RESPONSE 
El * El ♦ El 3 

X2 a El - B i*Xl - B0*X0 ♦ Ell 
El(I) IS SIGNAL AT A POINT IN FILTER 
E2U) IS SIGNAL AT A POINT IN FILTER 
IFIJJ.EQ.1) Eli 1 1 +2 ) =X2 
IF(JJ.EO.l) E2( II*2> = El 
EO IS FILTER OUTPUT 

EO * A2*EI-MAi-Bl*A2)*XH-I A0-B0*A2)*X0*EI2 
Ell* 0. 

E I 2 = 0. 

E 13 = 0. 

XO = XI 
XI » X2 

IFIJJ.EQ.I) GO TO 201 
I F I J J • EQ« 3 ) GO TO 202 
IFIJJ.EQ.4) GO TO 203 
F( J) IS X<N1~I)H< n 
El IS XIN1-I) 

YO IS HII) 

F I 3 ) a Y0*E 1 1 N1 + 3-I I ) 

F(A) = Y0*El(Ni«-2-II)*(-Bl) 

F < 5 ) a YO*E l ( N !♦ 1— 1 1 )*<— BO) 

F l 6 ) * F ( 4 ) ♦ F ( 5 ) 

EE II ) USEO TO CALCULATE RMS ERROR 
ERII) USED TO CALCULATE RMS ERROR 
EMI I) USED TO CALCULATE MAXIMUM ERROR 
DO 2 1=3,6 
EEin a ee i n FII) 

ERII) = ERII) ♦ FC I ) **2 
2 EMU) * EMI I ) M ABSIFI I ) ) 

GO TO 201 

202 F ( 7 ) * YO*E2(Nl+3-II )#A2 

F 1 8 ) a Y0*ElINl+2-in*IAl-ei*A2) 

F I 9 ) = Y0*E 1 (Nl+l-I I )*!A0-B0*A2) 

FI 10) = F I 7 ) ♦ F I 8 ) 

F I 2 ) = FC 9 > ♦ F 1 10 ) 

DO 3 1=7,10 
EE I I ) = EE(I) ♦ FID 
ERII) = ERII) + F I I ) **2 
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3 EMU ) * EMU) ♦: ABS{ F ( I ) ) 

EE ( 2 ) = EE ! 2 ) ♦: F ( 2 ) 

EH ( 2 ) = EK ( 2 ) + F { 2 ) **2 
EM( 2 ) = EM{ 2 ) * A8$< F ( 2 ) ) 

GO TO 201 

203 F ( 1 ) * YO*E2!Nl+3-II) 

EE ( 1 ) = EE ( i ) ♦. F(l> 

ERl 1) = ER! L) ♦ F ( 1 ) **2 
EMU) = EM ( 1) ♦; ABSCFU)) 

IFUl.LT.Nl) GO TO 201 

C SETS MULTIPLICATION ERROR TO ZERO IF COEFF ICIENT= 1 

IFIB1.E0.1.) EE!4)=0. 

IF1B1.EQ.1.) EM!4)*0. 

I F ! Bi • EU« 1. ) ER ! 4 ) =0» 

IF ! BO • EQ. 1 • ) EE ! 5 )=0# 

IF ( BO *E0* 1* ) EM ( 5 ) *0# 

IFtBO.EQ. 1. ) ER!5>=0. 

IFIA2.EQ.1.) EM ! 7 ) =0# 

IFCA2.EQ.I.) EE(7>*?0. 
lF!A2.fc0.l.) ER! 7 )=0» 

IF( ( Al-B1*A2) .EO. 1. ) EE!8)=0. 

IF!!Al-Bl*A2).EQ.l.) EM!8)*0. 

IF!(Al-Bl*A2).EQ.i.) ER!8)*0. 

IF((A0-B0*A2).6Q. 1.) EE!9)*0. 

IF!!A0-B0*A2).EQ.U> EM!9>*0. 

I F C I AO— BO* A 2 ) . EQ. 1 . ) ER!9>*0. 

00 50 1=1,10 

50 ERU) * E R ( I ) *XM2 ♦ EE U > **2*XM**2 
PRINT 300 

300 FORMAT! ’ - 1 » ’FORM’ *23X*’MAX ERROR 7X RM$ ERROR’ ) 

C CALCULATES ERRORS FOR INPUT QUANTIZER 

ERMI = SORT ( ER ( 1)*XM2> 

C CALCULATES ERRORS FOR OUTPUT QUANTIZER 

ERM2 = SORT ( ER ( 2 ) *XM2 ) 

PRINT 304, EM! U , ERMI 

304 FORMAT! •- • , • INPUT ’ • 15X , 2F 16. 5) 

PRINT 305* EM ( 2) , ERM2 

305 FORMAT! OUTPUT • , 14X,2F16.5) 

C SET 11*0 TO EXCLUDE INPUT ERROR POINT FROM TOTAL ERROR 

II = l 
11 = 0 

C SET 10=0 TO EXCLUDE OUTPUT ERROR POINT FROM TOTAL ERROR 

10 = 0 
10 = 1 

IFUUEQ.O) EE ! 1 ) =0 • 

IFU1-EQ.0) EM ( 1 ) =0 • 

IFUUEQ.O) ER I 1 ) =0 • 

IF(IO.EQ.O) EE ! 2 ) =0* 
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I F ( 10.60*0) EM|2>=0. 

IFU0.6Q.0I ER I 2) =0. 

EER IS USED IN CALCULATNG RMS ERROR FOR TRUNCATION 
EMM IS TOTAL MAXIMUM ERROR FOR FILTER FORM 
EER = 0. 

EMM = 0. 

ERR = 0. 

DO 10 1=1,10 
EMM = EMM+ EMI I ) 

10 ERR = ERR ♦ ER I I ) 

IFUT.NE. I) GO TO 40 
DO 30 1=1,9 
K * I + l 
DO 30 J=K ,10 

30 EER = EER ♦ EE ( !)*EEC J) 

EER = 2 • *EER 
C ERMS IS TOTAL RMS ERROR FOR FILTER FORM 
40 ERMS * SORT I ERR ♦ EER*XH*»2) 

PRINT 302, EMM, ERMS 

302 FORMAT! '-••■MODIFIED C ANON! CAL • , 2X , 2F 16 ♦ 5 > 

201 CONTINUE 
RETURN 
END 
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SUBROUTINE FIL2CEI|E0> 

CANONICAL FORM 

DIMENSION EH 150) ,E2C 1501 »£3l 150) *FU3) 
C0MM0N/C1/A2# Ai ,AO,B1,BO»XO»XI ,.X2»X3 t X4 
C0MM0N/C2/El,E2,E3,F t ER,EM,EE, I I , JJ,KK, 
X2 = El - B I*X 1 - BO*XO + til 


Ell * 0. 

IFC JJ.EQ. 1) E1UI+2) * X2 

EO = A2*X2 ♦! A l *X l + AO*XO ♦ EI2 

E 1 2 = 0. 

XO = XI 


XI a X2 

IFCJJ.EQ.l) GO TO 201 
IFCJJ.EQ.3) GO TO 202 
F C 3) = Y0*E 1 (N1+3-1 1 1 
F ( A) = Y0*E1CN1*2-I I)*C-B1> 
F ( 5 > = YO*El<NU-l-I I 1*(-B0J 
F ( 6 ) = F ( A ) ♦ FI 5 J 
DO 2 1 28 3* 6 

EE C II * EE ( I ) ♦: Fill 
ERU) = ERU) ♦ FU)**2 
2 EMU) * EMU) *. ABS(FUI) 

GO TO 201 

202 F ( 7) = Y0*El(Nl«-3-m*A2 
F { 8 ) a YO*E 1 INi*2~II)*Al 
F ( 9 ) = Y0*EHN1U-II)*A0 
F(10) = F (7 ) ♦ FC8) 

00 3 1=7, 10 


EEC I ) 
ER(I) 
EMC I ) 
IF (I I 
IFCBl 
IFCBl 
IFCBl 
IFCBO 
IFCBO 
IFCBO 
IFC A2 
IFCA2 
IFC A2 
IFCA1 
IFCAl 
IFCAi 
IFC AO 
IFCAO 
IFC AO 
00 50 


= EEC I- 1 
a ERU) 
= EMC I ) 
. L T . N 1 ), 

. EO* 1 • ) 
.EQ.l.) 
.EO. 1. ) 
.EO* 1. ) 

• EQ. I • ) 

• EQ.l.) 

• EQ* 1.) 
.EO. 1.) 

• EQ.l.) 

• EQ. 1* ) 

• EQ.l.) 
.EQ. 1.) 
.EQ.l.) 
.EQ. 1.) 
.EQ.l.) 

1 = 1,10 


+: fid 

* fU)**2 

♦ ABSC F ( I ) ) 
GO TO 201 

EE C 4 ) = 0. 
EMt 4 ) = 0. 
ERC4) = 0. 

EE ( 5 ) = 0. 
EMC5I = 0. 

ER C 5 ) = 0. 

EE < 7 ) = 0. 
EMC 7 ) = 0. 
ERt 7 ) = 0. 
EEC B I = 0. 
EMC8) = 0. 
ERC8) = 0. 
EEC9) = 0. 
EMC 9 ) = 0. 
ERC9I * 0. 


, ERC 13) ,EMCl3)tEEC 13) 
,X5,EIl,EI2,ei3 
N1,Y0, JT,XM,XM2 
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50 


ERCI) = ERU)*XM2 ♦ EE I I )**2*XM**2 
EER » 0. 

EMM = 0. 

ERR = 0. 

DO 10 1*1,10 

EMM = EMM ♦ EMt II 

10 ERR = ERR ♦ ER( I) 

I F C JT.NE. 1) GO TO AO 
00 30 1*1,9 
K = l+l 
DO 30 J=K ,10 

30 EER * EER ♦ EE ( I)*EE(J) 

EER * 2 •* EER 

40 ERMS * SORT (ERR ♦ EER*XM#*2) 

PRINT 302, EMM, ERMS 

302 FORMAT ( •CANONICAL*, 11X*2F16.5) 

201 CONTINUE 
RETURN 
END 
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SUBROUTINE FIl3IEItEO) 

C MODIFIEO DIRECT FORM _ 

DIMENSION Ell 150) *E2< 150)#E3U50)tFI 13) tER(13> *£M<13) *EE( 13> 
COMMON/ CI/A2, A 1 * AO, B 1 * BO * XO* X 1 *.X2 * X3,X4,X5,EI 1*EI2»,EI 3 
C0MM0N/C2/Ei,E2*E3fF»ER»EM*EEt II t JJ»KK,N1 *Y0, JT ,XM,.XM2 
X2 - A2*EI ♦ A i*X5 ♦ A0*X4 - X3+EI1 
Ell = 0. 

I F ( J J . EQ« 1 ) El< tI+2) = X2 
I F ( J J • EQ* 1) E2 ( 1 1 *2 ) =E I 
EG = X2 

X3 = 91*X2 * BO*X 1 
XI = X2 
X4 = X5 
X5 » El 

IF(JJ.EQ.l) GO TO 201 
F ( 3 ) = YQ*E 1 ( NI+3—I I ) 

F C A > = Y0*El(Ni + 2— II ) * I — Bl) 

F ( 5 ) = YO*EUNl-H-II)*<-BOa 
F ( 6 ) = FC4) + F I 5 ) 

F ( 7) = YO*E2(Nl+3-II)*A2 
F ( 8 ) = Y0*E2lNi+2-II)*Ai 
F I 9 ) = Y0*E2(Nl+l-in*A0 
F ( 10 ) = F ( 7 ) ♦ F ( 8 ) 

Fill) = F ( 9 ) ♦ F( 10) 

DO 2 1 = 3* 11 

EE< I ) = EE ( I ) ♦: FCI) 

ER I I ) = £R ( I ) 4? F(I)**2 
2 EM ( I ) = EM1I) * ABStFUl) 

IFUI.LT.NU GO TO 201 
IFI81.EQ.U) EE I A ) -0 • 

I F ( B1 *E0« 1 • ) EMI 4)=0# 

IFIBUEQ.l.) ERI4>*0. 

IFtBO.EQ. 1.) EE < 5 ) *0. 

IFIBO.EO.l.) EM(5)*0. 

IFIBO.EO.l.) ER(5)=0. 

I F I A2 *EQ. 1 • ) E0C7)=O. 

IFIA2.EQ.1.) EM ( 7 ) *0 . 

IF ( A2#£Q. 1* ) ER< 7 )=0# 

I F ( A 1 • EO* 1 « l £E< 8 )=0. 

IFIAWEQ.l.) £MI8)*0. 

IFCAi.EQ.l.) £R18>=0. 

IF (AO.EQ. 1* ) EEI 9 )=0« 

I F { AO. EQ. 1 • ) EMI 9 )=0. 

I F ( AO.EQ* 1*1 £RI9)=0. 

DO 50 1*1,10 

50 ER I I ) = ERII)*XM2 ♦ E E ( I ) **2*XM**2 
EER = 0. 

EMM = 0. 
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10 


ERR = 0. 

DO 10 1=1, 11 
EMM = EMM ♦ EMU) 

ERR = ERR ♦ ERill 
ZF< JT.NE. 1) GO TO 40 
DO 30 1=1,10 
K * I ♦ 1 
DO 30 J=K ,11 
30 EER = EER ♦ EE C l I *EE ( J ) 

EER = 2 •♦EER 

40 ERMS = SORTCERR ♦ EER*XM**2) 

PRINT 302, EMM, ERMS 

302 FORMAT (•-•,* MOO IF I ED D IREC T •, 5X, 2F 16* 5 I 
20! CONTINUE 
RETURN 
END 
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SUBROUTINE FlL4(El,E0) 

C MOOIFIED STANDARD FORM 

DIMENSION fcl (1501 ,E2( 150) ,E3< 150),FC13> ,ER( 13) ,EM(13) ,EE( 13) 

COMMON/ Cl/ A2, A l, AO, Bl, BO, XO * X l , X2 * X3 , X4 f X5 1 E I 1 , E I 2 * £ I 3 

COMMON/ C2/E 1,E2,E3»F,ER»EM,6E, II »JJ,KK,Nl f YO, JT, XM, XM2 

EO = A2*E I 4j XO 

X2 = AO*E I - B0*E0 ♦ Ell 

X3 = Al*fc I * XI - Bi*EO + E 1 2 

Ell = 0* 

E I 2 * 0. 

IFCJJ.EQ.l) £1(114-2} = X2 
I F < JJ.EQ. 1) E2( 114*2) = El 
IFIJJ.EQ.I) E 3 ( 1 1 4-2 1 = XI - Bl*EO 
XO * X3 
XI = X2 

IF( JJ.EQ. 1) GO TO 201 
IF< JJ.EQ. 3) GO TO 2 02 
F ( 3 ) * Y0*Ei(Nl+3-II ) 

F(4) = Y0*E2(Nl4-3-I I ) *A0 

F(5) * Y0*(Al*62(Nl4-2-in4-63(Nl4*2-II )4-A2*E2(NU3-U))*(-B0) 

DO 2 1=3,5 

EE ( I ) = EE ( I ) ♦ F ( I ) 

ER ( I ) = ER ( I ) 4= F( I )**2 

2 EMU) * EMU) ♦: ABSIFUM 
GO TO 201 

202 F ( 6 ) = Y0*E 3 ( N 14-3- I T ) 

F ( 7 ) = Y0*t2(Nl+3-I I >*Al 
F(8) = Y0*E2(Ni4-2-lI ) *A2 

F(9) = Y0*(E3(N14*1-II) 4- Al*E2(Nl4-l-H)*A2*E2(Nl4-2-im 
F ( 10 ) = YO*(Al*E2(Nl*2-in*E3(Nl*2-II)*A2*E2<Nl4>3-I 1 1 ) ( — S 1 > 
Ft 11) = F<6) ♦ F( 7) 

DO 3 1 = 6,11 

EE < I ) = EE ( I ) 4* F ( I ) 

ER(I ) * fcR ( I ) 4*. F ( I )**2 

3 EMU) = EMU) * A6S ( F ( I ) ) 

IFUI.LT.N1) GO TO 201 
IFtAO.EO. 1. ) EE(4>=0. 

I F ( AO • EO. 1.) EM<4)*0. 

I F ( AO • EQ. 1. ) ER < 4 ) =0 • 

IFiBO.EU.l.) EE ( 5 > *0. 

I F ( BO.EQ. 1. ) EM ( 5 ) =0 • 

IF(B0.E0.1.) ER( 5 ) =0. 

IF(Al.EQ.l.) £6 { 7 ) =0 • 

IF ( A 1 .EQ. 1. ) EM( 7 ) =0. 

IF( Al.EO. 1. ) ER ( 7 ) =0 • 

IFIA2.EQ.I.) EE ( 8 ) »0. 

IF(A2.EQ.1.> £M< 8 ) *0. 

IFCA2.EQ.1.) ER ( 8 ) =0 . 
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IFIB1.EQ. 1.) EE < 1 0 ) *0* 

I F < B 1 *£Q* 1 • ) EMI i 0 1 =0, 

IFIBi.EQ.l.) ER 1 103=0* 

DO 50 1=1, li 

50 ER C I ) = ER(I)*XM2 ♦ EE{I)**2*XM**2 
EER = 0. 

EMM * 0. 

ERR = 0. 

DO 10 1=1,11 

EMM = EMM 4- EMI I ) 

10 ERR = ERR ♦ ER ( II 

IF ( JT.NE. 1) GO TO 40 
DO 30 1=1,10 
K = !♦! 

DO 30 J = K » 1 1 

30 EER = EER ♦ EEI I)*EEIJ) 

EER = 2.*EER 

40 ERMS = SORT (ERR ♦ EER*XM**2) 

PRINT 302, EMM, ERMS 

302 FORMAT MODIFIED $ T ANDARD 1 *3X ,2F 16. 5) 
201 CONTINUE 
RETURN 
END 
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SUBROUTINE FIL5<EI*E0) 

C STANDARD FORM 

DIMENSION EII150) *E2U501 *E3C150)*FI 131 ,ERU3)fEMU3> *EE( 13) 

C0MM0N/C1/A2, Al f AO , B l , BO * XO, X l * X2, X3 * XA , X5, E I 1,EI2,EI3 

C0MM0N/C2/£i*£2*E3*F*ER*EM*6E*llfJJfKK*.Nl * YO* J T *XM*.XM2 

X2 = (A0-A2*B0-Bl*(Al-A2*Bi))*Et-Bl*Xl-B0*X3-*EIl 

XA = Xl+(Ai-A2«Bi)*EI+EI2 

EQ = X 3+A2*E I ♦£ 1 3 

Ell = 0, 

E 1 2 ^ 0. 

E I 3 = 0. 

IF(JJ.EO.l) El( 1 1 + 2 ) =X 2 
IF(JJ.EQ.l) E2I 1 1 ♦2) a E 1 
XI = X2 
X 3 = XA 

I F ( J J . EQ. 1) GO TO 201 
IF1JJ.EQ.3) GO TO 202 
IFUJ.EQ.A) GO TO 203 
F ( 3 ) = Y0*6 1 ( Nl+3- II) 

F ( A ) = Y0»El(Nl+2-II)*(-Bi) 

F I 5 ) = Y0*E2(Nl+3— II)*(AO-A2*BO— Bl*(Al-A2*8l)l 
F ( 6 ) * YO* ( ( A 1— A2 *B 1)*E2CN1 + 2~II )+El(Nl-fl — II) ) * I — BO ) 

FIT) = F I A ) ♦ F ( 6 1 

DO 2 1=3, 7 

EE ( I ) = EE ( I ) £ F( I) 

ER(I1 = ER ( I ) ♦ FU)**2 

2 EMU > * EMU ) * ABSIFU) 1 
GO TO 201 

202 F(8) = Y0*E2 (Nl+3— II)*CA 1-A2*B1 ) 

F ( 9 ) = F C 8 ) + Y0*E1 < Nl+2-If 1 

DO 3 1-8*9 

EE( I) = EE ( I) * F( I) 

ER C I ) = ER( I > «•' F ( I ) **2 

3 EMU) = EMU) ♦ ABSIFIO) 

GO TO 201 

203 F ( 10 ) = Y0*62(NU-3— I 1 )*A2 
EE< 10) = EE ( 10 ) ♦ F ( 10 ) 

ER( 10) = ERI10). ♦ FI 101**2 
EM ( 10) = EM MO ) ♦ ABS(F(10)) 

I F ( I I .LT.N1) GO TO 201 
IFIB1.EQ.1.) EE C A ) =0# 

I F ( B 1 *E0» 1* ) EM ( A ) *0 • 

IFIBI.EQ.I.) ER { A )=0 • 

IF U A0-A2*B0-B1*< Al-A2*Bl> J.EQ.l. ) EE<5)*0, 

IFC ( A0-A2*B0— Bl*< Al-A2*Bl> ).EQ.l. ) EM(5)=0. 

IF C C A0-A2*B0-6i*< A1-A2*B1) ) .ECUl. ) ER<5)=0. 

IFU Al-A2*B1 ).EQ. 1.) EE(8)*0. 

IF( ( A1-A2*B1).£Q. 1.) EM(8)*0. 
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IF! ! Ai-A2*Bl>.EQ.l*) Eft 1 8 ) *0* 
IF(BO.EQ.U) EE ! 6 ) =0* 

IF(BO.EQ.l.) EM ! 6 ) =0 • 

I F ( BO . fcQ* 1 • ) t R ( 6 ) ^0 • 

IF* A2.EQ. 1. ) EEU0>*0. 

IFCA2.E0.1.) EM!10>=0. 

IP(A2.EQ.l.) ER ! 10 I = 0 • 

00 50 1=1,10 

50 ERU) = ER(!)*XM2 ♦ E 6 ( I) **2*XM**2 

EER = 0. 

EMM = 0. 

ERR « 0. 

DO 10 1=1,10 

EMM = EMM «■ 6M( I) 

10 ERR = ERR ♦ ER ( I ) 

IF ( JT.NE. U GO TO AO 
DO 30 1=1,9 
K = 1 + 1 
DO 30 J=K ,10 

30 EER = EER ♦ EE (I ) *EE ! J ) 

EER = 2. ♦EER 

AO ERMS = SORT ( ERR ♦ EER^XM^2) 

PRINT 302, EMM, ERMS 

302 FORMAT! •-• , •STANDARD' , 1 2X , 2F 16* 5) 
201 CONTINUE 
RETURN 
END 
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SUBROUTINE FIL6IEI.E0) 

C PARALLEL FORM 

DIMENSION ELI 150) ,E2I 150) *E3< 150) ,F I 13) ,ER( 13) ,EM(i3) ,EE( 13) 
COMMON /C1/A2,A1,A0,B It BO, XO , X l , X2 , X3 , XA , X5 , £ I l,EI2,EI3 
C0MM0N/C2/E 1 ,E2,E3»F»ER,EMrfEE» 1 1 , J J,KK, Nl , YO, JT, XM, XM2 
C0MM0N/C3/Pi,P2»P3,P4,Gl,G2,G3,GA 
XI = G1*E I - P1*X0 ♦ Ell 

X3 = G2*EI - P2*X2 ♦ E 1 2 

EO = A2 *E I F XO ♦ X2 ♦ EI3 

Ell = 0. 

E 1 2 = 0, 

E 1 3 = 0. 

IF < JJ .EQ. 11 El ( I I +2) * XI 
IFIJJ.EQ. 1 ) E2I 1 1 ♦Z ) = X3 
IFIJJ.EQ. 1) E 3 ( 1 1 +2 I = El 
XO = XI 
X2 = X 3 

IFIJJ.EQ. 1) GO TO 201 
IFIJJ.EQ. 31 GO TO 202 
IFIJJ.EQ. A) GO TO 203 
F < 3 I = YO*E3(NL + 3-II)*Gl 
F I A) * Y0*El(Nl + 3-II ) 

F I 5 ) = YQ*E 1 ( N 1 ♦2— II 1*1— PIT 
DO 2 1=3,5 
EE I I ) = EE C I ) F FID 
ERI I) = ERI I ) F FI I)**2 

2 EMU) = EMII) *■ ABSI F (1)1 
GO TO 201 

202 F I 6) = YQ*E3 I Nl + 3— 1 1 ) *G2 
F I 7) = Y0*E2I Nl + 3- II) 

F I 8) = YO*E2(Nl+2-in*(-P2> 

DO 3 1=6,8 

EE I I ) = EE ( I ) ♦ FII) 

ERII) = ERII ) F FI I )**2 

3 EMII) = EMI I) ♦ A0SIFI I ) ) 

GO TO 201 

203 F ( 9 ) * YO*(El(im2-Vn+E2<Nl+2-IIM 
F I 10 ) = Y0*E3I N1F3-I I )*A2 

00 8 1=9, 10 

EE I I ) = EEI l) ♦ FII) 

ERII) = ERII) F FII)**2 
8 EMII) * EMII) F ABSI F I I I ) 

IF! II.LT.N1) GO TO 201 
IFIGl.EQ.i.I EE I 3 ) =0 . 

IFIG1.EQ.1.) EMI 3 > =0. 

IFIGl.EQ.i.I ER I 3 )*0. 

IF(Pl.EQ.U) EE I 5 ) *0. 

IFIP1.E0.1.) EM ( 5 1=0 • 


61 



IF(Pl.fcQ.U) £R ( 5 } =0 • 
l F ( G2 • E Q* i • ) EE ( 6 ) *0 • 

IHG2.EQ. 1. ) EM ( 6 ) =0 • 

I F ( G2 *EQ* 1* ) ER(6>*0. 

IFIP2.EQ.W) EE ( 8 ) =0. 

IF(P2.EQ.l.) EM<8)*0. 

I F < P 2 • E0« 1 • ) ER ( 8 ) = 0 • 

IFIA2.EQ.1.) E £ ( 1 0 > = 0 . 

I F ( A2 • EQ • l . ) EMdOMO. 

IHA2.EQ.I.) Eft ( 1 0 ) -0 • 

DO 50 1=1*11 

50 Eft ( I ) = ER(U*XM2 + EE (I )**2*XM**2 


EER = 0, 

EMM = 0. 

ERR = 0. 

DO 10 1*1,10 

EMM = EMM + EM(I) 

10 ERR * ERR ♦ ER( I) 

IF( JT.NE. 1) GO TO 40 
DO 30 1*1,9 
K = IU 
DO 30 J=K ,10 

30 EER * EER ♦ EEC I)*EE( J) 

EER * 2 ®*EER 

40 ERMS * SORT < ERR + EER*XM**2) 

PRINT 302 , EMM, 6RMS 

302 FORMAT I 8 -« ,« PARALLEL • , 12X,2F16.5> 

201 CONTINUE 
RETURN 
END 
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SUBROUTINE FIL7(EI,EO) 

C XI FORM 

DIMENSION El 1150) ,E2I150) ,E3<150) ,F(13) ,ER(13) ,EMU3) ,EE( 13) 
C0MM0N/Ci/A2tAl # A0fBl,B0,X0,Xl f X2,X3,X4,X5, £1 1,EI2,EI3 
COMMON/ C 2/ El,E2»E3,F,ER»EM^EE,II»JJ*KK,Nl,YG*JTtXM* XM2 
C0MM0N/C3/P l»P2,P3,P4*Gi*G2,G3*.G4 
XI = 2.*G3*6I - P3*X0 - P4*X2 4 Ell 
X 3 = 2.*G4*6I - P3*X2 ♦ P4*X0 4 612 
EO = X2 + A2*E I 4 E 1 3 
Ell = 0. 

£12 * 0. 

£13 * 0. 

IF< JJ.EQ. 1) Ell II + 2) = XI 
IF( JJ.EQ. U E 2 ( 1 1 +2 ) = X3 
IFCJJ.EQ.i) 63(1142) * El 
XO = XI 
X 2 = X 3 

IF< JJ.EQ. 1) GO TO 201 
IF (JJ.EQ* 3) GO TO 202 
I F ( J J • E Q« 4 ) GO TO 203 
F ( 3 ) = YO*E l(Nl+3-II > 

F ( 4 ) = YQ*E1(N142-II ) *<-P3 ) 

F ( 5 ) = Y0*E3(Ni43-II)*2.*G3 
F ( 6 ) = Y0*E2(N142-U)*(-P4) 

F ( 7 ) = F( 5 I 4 F ( 6 ) 

DO 2 1 = 3,7 

EE ( I ) = EE ( I) 4 FUI 
ERU) = ERI I) 4; F(I)**2 

2 EMU) = EMU) 4 ABS(FU)) 

GO TO 201 

202 F ( 8 ) = Y0*E2(N143-II) 

F{9) = Y0*E2(Ni42-in*t-P3) 

F ( 10 ) = Y0*E3CN143-I I)*2.*G4 
Fill) = Y0*8l I N 14.2*1 I ) *P4 
FI12) * F ( 10) 4 F ( 1 1 ) 

00 3 1=8, 12 

EEII) = EE II ) 4 FID 

ERU) = ERU) 4 . F(l>**2 

3 EMU) = EMU) 4 ABS(FU)) 

GO TO 201 

203 FUJI = YO*E3|Nl + 3-II ) *A2 
EE I 13 ) = EE M3). 4 F( 13) 

ERU3) = ERI 13). ♦ FI 13)4*2 
EM I 1 3 ) = EMI 13) 4 AQS ( F( 13) ) 

I F ( 1 1 .LT.N1) GO TO 201 
IFIP3.EQ.1.) EE 1 4 ) = 0 • 

I F ( P3 • LO. 1 • ) EM|4)=0. 

IFIP3.EQ.1.) ER ( 4 ) *0 » 
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I M < 2 • *G3 ) * 60. 1 • ) EE(5)=0* 

I F ( { 2 «*G3 ) . 60 - I . ) EM{5)=0. 

IF( I 2**G3 I . EQ. 1* ) ER I 5 ) =0. 
IFIP4.EQ.1.) EE ( 6 ) =0# 

If(P4*EQ.l.) EM ( 6 ) =0 • 

IF(P4.EQ.1.) ERt6>=0. 

IF { P3 • EQ* 1 * ) EE ( 9 )=0* 

I F ( P3 • EQ« 1. ) EMI 9 ) =0 • 

IFIP3.EQ. U ) ER ( 9 I =0 • 

IF ( ( 2 **G4 ) • EC « 1 . ) EE(10)=0. 
IF(<2.*G4).EQ.K) EM ( 1 0)=0* 
IFC(2.*G4).e0.1.) ER I 10)=0. 
IFIP4.EQ.1.) EE(ll)*0. 

I F ( P4 *EQ« 1 • ) EM { 1 1 ) = 0 • 

IF(P4.EQ. 1. ). ER< 1 15=0. 

IFIA2.EQ.1.) EE ( 1 3 ) =0 • 

I F ( A2 • EQ. 1 . ) EM(13)»0. 

IF( A2 «EQ» l. ) ER ( 1 3 } = 0 . 

DO 50 1=1,13 

50 ER(I) = ER ( L ) *XM2 ♦ EE < I) **2*XM**2 
EER = 0, 

EMM = Oo 
ERR = 0. 

DO 10 1=1,13 
EMM = EMM + EMIII 
10 ERR = ERR ♦ ER( I) 

IFIJT.NE.il GO TO 40 
DO 30 1=1,12 
K = 1*1 
DO 30 J=K , 13 

30 EER = EER ♦ EE( I)*EE( J) 

EER = 2 • *EER 

40 ERMS = SORT ( ERR ♦ EER*XM**2> 

PRINT 302, EMM, ERMS 
302 FORMAT! • XI • , 18X,2F 16. 5) 

201 CONTINUE 
RETURN 
END 
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SUBROUTINE FIl8<EI,EO) 

XII F0R ^ . , 

DIMENSION EIU50) .E2U50) »E3(150) *F(13> *ER( 13) t EMU3> ,EE( 13) 

C0HM0N/Cl/A2»Alt A0*Bl»B0*X0*Xl,X2,X3,X4 f X5*EI i*EI2,EI3 

C0MM0N/C2/EI,E2,E3, E,ER,EM#EE, 1 1 , JJfKK^Ni f Y0* JT»XM t XM2 

C0MM0N/C3/P 1* P2* P3(>P4*Gl* G2*G3*G4 

XI = El - P3*X0 ♦ P4*X2 + Eli 

X3 - -P3*X2 - P4*X0 ♦ E 1 2 

60 = 2.*G4*X0 - 2.*G3*X2 ♦ A2*EI + EI3 

Ell « 0. 

E I 2 * 0. 

E I 3 = 0. 

IP(JJ.EQ.I) EH 1 1 +2 ) = XI 
IF< JJ.EQ. 1) E2< 1 1 +2 ) = X3 
IF< JJ.EQ. 1) E3C 1 1 +2) = El 
XO = XI 
X2 * X 3 

IF( JJ.EQ. 1) GO TO 201 
IF (JJ.EQ. 3) GO TO 202 
IFCJJ.EQ.4) GO TO 203 
F ( 3 ) = Y0*EKNU3-II) 

F ( 4 ) = YO*El(Nl+2-II)*(-P3) 

F ( 5 ) * Y0*(63(Nl+3-II )+P4*E2(Nl+2-II) ) 

F ( 6 ) = Y0*E2 ( N1+2—I I )*P4 
DO 2 1-3*6 

EE ( 1 ) * EE ( I ) ♦ FID 
ER ( I) = ER ( I ) ♦ F(I)**2 

2 EMU) = EMU) ♦: A BS ( F ( I ) ) 

GO TO 201 

202 F ( 7 ) = Y0*E2(Ni + 3-I I) 

F I 8 ) = Y0*fc2(Nl + 2-II )*(-P3) 

F ( 9 ) = Y0*E1(NI*2-II )*P4 

DO 3 1=7*9 

ECU) = Efc(I) ♦' FII) 

ERU ) = 6R( I) FU»**2 

3 EMU) = EMU) ♦ ABSlFCin 
GO TO 201 

203 F ( 10 ) = Y0*E2(NH2-II )*2.*G3 
Fill) » Y0*El(Nl+2-I I )*2.*G4 
F ( 12 ) = F ( 10 ) ♦ Fill) 

F ( 13 ) = YO*E3(NH3-II )*A2 

DO 8 1=10*13 

EE < I ) = EEC I- ) ♦' FU) 

EK( I ) = ERU) +i FC I)**2 
8 EMU) = EMU) ♦ ABSCFCIM 
IFCII.LT. Nl). GO TO 201 
IFCP3.EQ.1.) E&( 4 ) =0 • 

IFCP3.EQ.I.) EM ( 4 ) *0 • 
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IFCP3.EQ.1.) ERC4|=0. 

IFCP4.EQ.1.) 6B( 6 )=0« 

IFIP4.EQ.U1 EMC 6 ) =0* 

IFCP4.EQ.I.) ER:(6>*0. 

IFCP3.EQ.1.) EE l 8 ) =0. 

I F C P 3 • EQ. 1.5 EMC 8 )*Q. 

IE ( P3 . EQ. 1 • ) ERC 8 )=0. 

IFCP4.EQ.1.) EM ( 9 ) =0 • 

I F ( P4 . EQ. 1 • ) EEC 9)=0. 

IF C P4.EQ. U ). HR C 9 ) = 0 • 

IFC C2.4G3U6Q.1.) EEC 10 )=0. 

IF ( ( 2 .*G3 ) • EQ . 1 . ) EMC i 0 ) *0. 

IFC C2.*G3).BQ.l.) ERC10}=0. 

I F C ( 2 .*£4 U6Q. 1. ) EEC 1 1 ) =01 
IFC (2.*G4)«£Q.lo) EMC 1 1 5 = 0# 

IF( (2.*G4) .EQ.l.) ER(ll)*0. 
IFCA2.6Q.1.1 EEC 13) =0. 

IFCA2.EQ.1.) EMC 13)*0. 

I F C A2 • EQ. 1 . ) ER ( I 3) *0 . 

00 50 1=1.13 

50 ERCI) * ERCI) *XM2 ♦ EE C I ) **2*XM**2 

EER = 0. 

EMM * 0. 

ERR = 0. 

DO 10 1=1.13 

EMM = EMM + EMC I ) 

10 ERR = ERR ♦ ERCI) 

IFCJT.NE. I) GO TO 40 
00 30 1=1,12 
K = in 
00 30 J=K ,13 

30 EER = EER ♦ EEC I )*EEt J) 

EER » 2.*EER 

40 ERMS * SORT C ERR ♦ EER*XM**2) 

PRINT 302, EMM, ERMS 
302 FORMAT* 17X.2F16. 6) 

201 CONTINUE 
RETURN 
ENO 
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